Orientational Ordering and Dynamics of Rodlike Polyelectrolytes 



Hossein Fazli,^ Ramin Golestanian,^ and Mohammad R. Kolahchi^ 

''institute for Advanced Studies in Basic Sciences, Zanjan 45195-159, Iran 
(Dated: February 6, 2008) 

The interplay between electrostatic interactions and orientational correlations is studied for a 
model system of charged rods positioned on a chain, using Monte Carlo simulation techniques. It is 
shown that the coupling brings about the notion of electrostatic frustration, which in turn results in: 
(i) a rich variety of novel orientational orderings such as chiral phases, and (ii) an inherently slow 
dynamics characterized by stretched-exponential behavior in the relaxation functions of the system. 
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I. INTRODUCTION 

Solutions of highly charged polymers are known to de- 
velop novel structural properties due to the interplay be- 
tween electrostatic interactions and entropic effects 
For example, it has been shown that for most charged 
biopolymers such as DNA, filamentous (F-)actin, and 
various viruses Q , electrostatic correlations in the vicin- 
ity of the macroions caused by multivalent counterions — 
ions of opposite charge — can lead to like-charged attrac- 
tion 3 . This attraction most often destabilizes the poly- 
electrolyte solution and leads to the formation of col- 
lapsed bundles. 

In a recent experiment with F-actin, Wong et al. ob- 
served that when the density of multivalent counterions 
is not yet sufficient to trigger complete collapse, the uni- 
form solution can become unstable and phase separate 
into coexisting low- and high-density phases [4] . This ex- 
periment is perhaps a most direct manifestation of the 
peculiar behavior of polyelectrolytes in high density so- 
lutions: the dense phase is characterized by multi-axial 
liquid crystalline behavior 0] , as well as exceedingly slow 
dynamics that leads to the formation of an actin gel [3]. 
A similar self-assembly pattern has been observed in the 
structure of the nuclear lamina — a charged filamentous 
integral membrane protein network that provides a cy- 
toskeletal support for the nuclear membrane 

Another interesting example is the observation of slow 
modes in the dynamics of various high-density polyelec- 
trolyte solutions in the absence of salt, where it is sug- 
gested that correlations cause the formation of starlike 
complexes with considerably slow dynamics 0, [l] . These 
examples raise the issue that strong many-body Coulom- 
bic correlations for extended (non-pointlike) charged ob- 
jects could lead to novel features that are poorly under- 
stood. 

Based on such experimental evidence, one can broadly 
identify two characteristic features in this class of prob- 
lems: the novel orientational ordering, and the slow dy- 
namics. As a first step towards understanding such prop- 
erties, one should note the frustration that is inherent 
structurally in the collection of like-charged rodlike ob- 
jects: As they come together, they reorient themselves 
to minimize the total energy of the system. Here, the 
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FIG. 1: The schematic configuration of charged rods. The 
centers of the rods of equal length L are fixed on a chain of 
spacing a. Each rod has two rotational degrees of freedom, 
and a uniform distribution of point charges. 



long range nature of the interaction presents them with 
a multitude of configurations that are good minimum en- 
ergy candidates and yet are far apart in the configuration 
space. This picture is partly supported by the fact that 
there are so many of them. This picture is supported by 
the fact that in the presence of (a considerable amount 
of) salt the above effects are washed out, as screening can 
eliminate the frustration by reducing the effective range 
of the mutual repulsion. 

Here, we set out to study the collective behavior of 
polyelectrolytes at close separations using Monte Carlo 
(MC) simulation, with the primary goal of understand- 
ing the complexity that electrostatic frustration can bring 
about. We consider a most simple model system of simi- 
larly charged rigid rods of equal length L with their cen- 
ters forming a regular lattice of spacing a on a chain, as 
shown in Fig.[T] We find a variety of collective behaviors, 
and map out the phase diagram for the system as a func- 
tion of temperature and lattice spacing (inverse density of 
the rods), as shown in Fig. O In the low-density regime, 
the rods are found to order in a staggered way upon de- 
creasing the temperature, such that they are perpendicu- 
lar both to their neighbors and to the axis that connects 
their centers, through a two-stage crossover transition. 

For intermediate densities, we find that a low- 
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FIG. 2: Phase diagram of the system in the a — T plane. 
Dashed hne: The crossover transition hne from a disordered 
regime to partially ordered one [region (a)]. Dash-dotted 
line: The crossover transition line from a partially ordered 
regime to the staggered ground-state of the low-density regime 
[region (b)]. Solid line: First order transition line from a 
partially ordered regime to the chiral ordered phase of the 
intermediate-density case [region (c)]. Hashed region: The 
domain in which a multitude of transitions occur to phases 
with various kinds of orderings. Region (d) : The "sea urchin" 
phase. Schematic configuration of rods corresponding to the 
different parts of the phase diagram is indicated in the lower 
panel. The values of the various transition points indicated 
in the phase diagram are = 0.401/, ai = 0.10 L, and 
o* ~ 0.02 L. 



in Sec. IIVI Some discussions and concluding remarks are 
presented in Sec. |Vl Finally, a variant of the model that 
corresponds to polyelectrolyte combs is discussed in the 
Appendix. 



II. THE MODEL 

We use the Metropolis algorithm to simulate a system 
of Nr rods with their centers fixed on a chain, and only 
their rotational degrees of freedom {9, (p) left to explore 
. We assume that Nq (odd) point-charges of the same 
sign are distributed symmetrically on each rod (Fig. [1]). 
The screened Coulomb interaction energy of the system 
of rods can be written as: 
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where e is the dielectric constant of the medium, q — 
Q/Nq where Q is the overall charge of the rod, rmi(r„j) 
is the position of charge z(j) on rod m{n), and is the 
Debye screening length. Note that the energy expression 
is invariant under the local transformation — > tt — and 
(/) — > (/)-f TT for each rod. We use free boundary conditions 
to avoid complexities arising from the Ewald-summation 
of extended charge distributions [IC^]. 

For each value of the lattice constant a, we start the 
system from a random configuration of the rods at a 
high temperature (where the system is in the disordered 
phase) and gradually cool down toward lower tempera- 
tures. At a given temperature, we start the calculation 
of various thermal averages after equilibration and take 
the last configuration at each temperature as the initial 
configuration for a slightly lower temperature. We have 
run a wide variety of simulations for different values of 
Nr, Nq, and a, and found that all the qualitative results 
are robust. 



temperature ordered phase appears through a first-order 
phase transition, in which neighboring rods have a twist 
angle in addition to the ninety degrees of the low-density 
phase. In the high-density regime, a hierarchy of differ- 
ent phases is observed depending on the lattice spacing, 
where there is a multitude of twist angles, a periodic 
structure with a basis, and out-of-plane arrangements 
of the rods. We also study the equilibrium relaxation 
properties of the system and find that the appropriate 
order parameter in the system has an anomalous relax- 
ation characterized by a stretched-exponential behavior, 
which is a feature also seen in systems with frustration. 

The rest of the paper is organized as follows. In Sec. 
im we introduce the model and outline the simulation 
technique used. Section Hill is devoted to the description 
of various parts of the phase diagram, which is followed 
by discussions on the relaxation dynamics of the system 



III. THE PHASE DIAGRAM 

To capture the essence of the electrostatic frustration, 
we first consider the salt-free case that corresponds to 
K — 0, and comment on the effect of screening later. In 
Fig. [51 we show the phase diagram of the system in the 
plane of temperature and lattice constant for Nr = 100 
and Nq — 3. 

In the low-density regime, where a is greater than a 
critical value Oc = 0.40 L, the rods appear to order on de- 
creasing the temperature through a two-stage crossover 
transition. The first crossover (denoted by the dashed 
line in Fig. ^ takes the system to a partially ordered 
phase in which the angular degree of freedom, 9, freezes 
to ^ , and the rods fluctuate freely only within the result- 
ing parallel plates [region (a) in Fig. [5]. The ordering in 
the 4> degree of freedom occurs through another crossover 
at a lower temperature (denoted by the dash-dotted line 
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FIG. 3: The specific heat as a function of temperature T [in 
units of q^/{ek-BL)\ for (a) a = L and (b) a — 0.351/, with 
three different sizes iV^ = 50, 100, 200 and iV, = 3. 

in Fig. [2]). The crossover transition temperatures are de- 
termined from the specific heat of the system as a func- 
tion of temperature, which develops two rounded peaks. 
In Fig. [3K, the specific heat of the system is plotted 
as a function of the temperature T, which is made di- 
mensionless by the combination q^/[£k^L), at a = i 
for three different sizes Nr = 50, 100, 200 and Nq = 3. 
Each data point is the average over five independent runs 
starting from different initial random configurations. For 
each temperature, 1.5 x 10^ MCS are used with the first 
7.5 X 10^ MCS being for equilibration. As can be seen 
from Fig. [3^, the height of the peaks in the specific heat 
do not change appreciably for various values of the sys- 
tem size Nr = 50, 100, 200. 

The ground-state configuration corresponding to this 
regime is given as [see Fig. Od] : 

Om = and = m X -, (2) 

where m indicates the position of the centers of rods on 
the z axis (zm = ma) . This is in agreement with the two- 
body minimum energy configuration of the rods [11]. A 
multipole expansion of the electrostatic energy expression 
in Eq. ([T|), which is presumably a good approximation 
in the limit of a L, reveals that the system of rods 
effectively behaves as a collection of quadrupoles in an 
external staggered electric field; hence the perpendicular 
ordering described by Eq. This picture also agrees 
well with the lack of genuine phase transitions in this 
regime, as it is expected from quadrupoles in a symmetry 
breaking external field. 

In the intermediate-density regime, corresponding to 
a* < a < Uc where a* ~ 0.025 L, lowering the tem- 
perature for each value of the lattice constant causes a 
crossover transition from disordered to partially ordered 
regime, which is followed by a phase transition (solid line 
in the phase diagram in Fig. [2]) to low-temperature phases 
having chiral order. A plot for the specific heat in this 
regime, corresponding to a = 0.35 L, is shown in Fig. ^jp 
for system sizes = 50, 100, 200 and Ng = 3. The data 
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FIG. 4: The thermal average of energy per site [in units 
of q'^/{eL)] as a function of the temperature [in units of 
/(ekBL)] upon successive cooling and heating of the sys- 
tem, for a = 0.35 L, A^^ = 100, Nq = 3. The self-energy of 
the staggered ground-state configuration E+ is subtracted off 
from the energy. The hysteresis loop can be observed, which 
indicates a first order transition from the partially ordered 
regime to the chiral phase. 

points are the averages of five independent runs starting 
from different initial random configurations. For each 
temperature 2 x 10^ MCS are used, with 10^ MCS being 
for equilibration. It can again be seen that the specific 
heat is not sensitive to the size of the system. 

To understand the nature of the phase transition, we 
examined the energy of the system upon successive cool- 
ing and heating. In Fig. [H the thermal average of energy 
per site as a function of the temperature [that is made 
dimensionless using q- /{ek-BL)] is plotted for a = 0.35 L 
corresponding to successive cooling and heating of the 
system. In this curve, 7.5 x 10^ MCS are used for equili- 
bration and 7.5 x lO'^ MCS for ensemble averages, and it 
corresponds to Nr = 100, Nq — 3. As can be seen from 
Fig. m the average energy of the system shows hysteresis 
around the transition temperature, which is a signature 
of a first order phase transition. One can also define 
an order parameter that could distinguish between the 
two phases and thus monitor the phase transition as a 
function of temperature. A suitable choice for the order 
parameter could be 

n— 1 

where A„(i) = (^„+i(i) — (^„(t) — ^. Figure [5] shows a plot 
of this order parameter versus temperature [in units of 
q'^ /{ek-BL)] for a system of Nr — 200 and A^g = 3 at a = 
0.35 L. In this plot, 10^ MCS are used for equilibration of 
the system and 10^ MCS for thermal averages. It can be 
seen that the order parameter develops a discontinuous 
jump across the transition temperature. This is another 
indication that the transition across the solid line in the 
phase diagram of Fig.[2]is first order. These observations, 
however, should be complemented with more systematic 



4 




FIG. 5: The order parameter M [defined in Eq. ((Sjl] versus 
temperature [in units of / {eksL)] for a system of Nr = 200, 
iV, = 3 at a = 0.35 L. 

investigations using energy histogram methods [H, [l^l . 

Below this transition, the system of rods appears to 
develop a rich variety of orderings. For ai < a < 
where ai = 0.1 L [region (c) in the phase diagram in 
Fig. [2] , the array of rods orders in a chiral phase with 
the corresponding ground-state configuration given as: 

6l,„ = ^ and m x + Ag(a)y (4) 

Since the twist angle A is small in the vicinity of the tran- 
sition line, we can simplify our description by assuming 
that A itself serves as the corresponding order parameter 
for the phase transition from (a) and (b) to (c). While 
its value is observed to undergo a finite jump just below 
the transition line, the ground-state value Ag(a) is found 
to vanish continuously at Cc, implying that the termi- 
nation point of the transition line at zero temperature 
corresponding to a = Ac is a critical point. 

For a more direct study of the transition in the ground- 
state configuration at ac, we calculate the exact energy 
per site for the system using the energy expression given 
in Eq. ([T]) as a function of A for the configuration given 
in Eq. This means that we freeze the charged rods 
into the configurations defined by Eq. (g]) and calculated 
in a point by point double summation the electrostatic 
energy of the configuration for each value of A. This 
energy function, which is plotted in Fig. [H] for different 
values of a, admits a Landau form and the correspond- 
ing minimum energy solution Ag(a) vanishes at ac as 
Ag = |ac - a\l^ with 13 = 0.500 (see inset of Fig. We 
have checked that the value of (i does not depend on Nq 
(even up to Nq = 1001), confirming the expectation from 
Landau theory that /3 = 1/2. 

In the intermediate-density regime when the rods are 
brought closer to each other, different configurations that 
are close in energy appear. This changes the behavior of 
the system at ai where the curve of energy versus A 
starts to have four degenerate minima. For a* < a < ai 
with a* ~ 0.025 L (hashed region in Fig. [2]), while the 



FIG. 6: A plot of the ground-state energy per site as a func- 
tion of A for a = 0.3 L, a — ac — OAL, and a — 0.6 L, with 
Nr — 40 and Nq — 3. The value of A which minimizes the 
energy of the system vanishes continuously at a — a^. In- 
set: The log-log plot of Ag versus (oc — a) gives a slope of 
l3 = 0.500. 



rods still respect the planar ordering, a multitude of 
transitions occur taking the system through a variety of 
minimum energy configurations. For example, upon in- 
creasing the density, a phase appears in which successive 
cross-like structures formed by two neighboring rods are 
rotated relative to each other around the z-axis by an 
angle A', which corresponds to a chiral ordering with a 
basis. Finally, in the high-density region where a < a*, 
the rods are no longer confined in the planes. While 
in this regime one also finds various kinds of orderings, 
eventually the rods tend to approach a spherical configu- 
ration that resembles a "sea urchin" [see Fig.[2I^d)], as a 
approaches its smallest possible value oq that is set by the 
thickness of the rods. Note that only a very small choice 
of the cutoff length, such as our choice of oq = 10~^ L, 
can allow for these complicated high-density phases to ex- 
ist. Moreover, unlike the lower density cases, this part of 
the phase diagram depends very sensitively on the choice 
of charge distribution on the rods, namely the value of Nq 
as well as whether the charges are smeared or discrete. 



IV. RELAXATION DYNAMICS 

The low-temperature chiral ordering (in either of its 
various forms) is not easily achieved in its complete form 
when the system is cooled down from the disordered 
phase to temperatures below the transition line. This 
means that the equilibrium dynamics of the system in- 
volves a rather slow process, which we try to identify and 
study here. 

For ai < a < ac, the array of rods is decomposed into 
several domains of the chiral structure with both positive 
and negative values of A, corresponding to the degener- 
ate minima of energy shown in Fig. [51 The neighboring 
left-handed and right-handed domains are regions that 
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FIG. 7: Equilibrium relaxation function of the system for 
a = 0.25 L and various temperatures, reported in units of 
(f / {ekB,L). The nonlinear fits to the long-time behavior of 
the curves are of the stretched-exponential form. Inset: A 
plot of 7 as a function of temperature. 



typically contain just a few rods, and are linked via what 
can be thought of as kinks. By definition, the kinks anni- 
hilate when they meet each other or when they reach the 
boundaries of the system. We observe that the density 
of the kinks is more or less fixed to about one kink for 
every 40 rods or so, irrespective of the size of the system. 

These kinks are observed to have very slow dynamics, 
presumably due to the electrostatic frustration, so that a 
considerably long MC time is needed before all the kinks 
are annihilated and one of the domains spans the entire 
system. The MC time needed for the annihilation of the 
kinks increases with increasing the system size. We also 
note that when the energy has four minima (for a just 
below ai) and two kinds of chiral structures are possible, 
different kinds of kinks are observed depending on the 
type of the two neighboring domains that are linked. 

To put the study of the slow dynamics in a more 
quantitative framework, we probe the equilibrium relax- 
ation properties of the system by measuring the auto- 
correlation function for A„(t), which is defined as: 



C{t) 



(5) 



To obtain C(t) at each temperature T, we start the sys- 
tem from a random configuration at high temperatures 
and quench it to T. After running the system for a wait- 
ing time, tw, we use the system configuration 0„(i) in 
all time steps from tw to tw +tav to calculate C'(t) using 
Eq. (O. For long enough waiting times {t^ ~ 2 x 10^ 
MC steps for Nr — 100)), the relaxation function C{t) is 
independent of tw as well as the system initial conditions, 
indicating that there is no difference between averaging 
over time [Eq. ([S])] and averaging over different initial 
conditions of the system. Note that the waiting time is 
typically the time needed for the kinks to annihilate. 

In Fig. [3 we show the system relaxation function with 
logarithmic time scale for different values of tempera- 
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FIG. 8: First order transition line of the phase diagram for 
different values of the screening parameter. 
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q — 3. Using nonlinear fits with three free parameters, 
we find that the long-time part of the relaxation function 
is well represented by the stretched-exponential form, 
Co exp[— (t/r)^]. The exponent 7 saturates to a constant 
of the order of 0.9 in the high temperature regime, while 
it drops at lower temperatures where the relaxation is 
exceedingly slower. The stretched-exponential form for 
the relaxation is a feature that is often encountered in 
the dynamics of frustrated systems 
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V. DISCUSSION AND CONCLUSION 

In a most simple realization of the system, we have 
shown that charged linear objects in close separation can 
develop novel collective orientational correlations as well 
as an unusually slow dynamics, due to the frustration in 
the orientational degrees of freedom. While in our simple 
model system we assume that the rods are positioned in a 
lattice, we expect that removing this constraint will cause 
the orientational correlations to manifest themselves as 
liquid crystalline structures such as those observed in Ref. 
[J]. Geometrical frustration induced by the dimension- 
ality of the arrangement could also affect the orienta- 
tional ordering, and it remains to be seen what kind of 
arrangement is favored by the liquid system itself, i.e., 
what is the coupling between the positional- and the ori- 
entational correlations in this frustrated system. 

Let us examine the effect of screening by added-salt in 
the above results. For a low concentration of added salt, 
the screening is weak and we expect that the results re- 
ported here will hold true so long as the screening length 
is the largest length in the system. As the salt concen- 
tration is increased, however, the Debye length begins 
to compete with the other length scales and this could 
change the phase behavior of the system. 

We have performed simulations with the screened 
Coulomb interaction of Eq. ([1]) using different values 
for the Debye screening length. Wc find that the general 
structure of the phase diagram of Fig. [2] persists, while 
the size of the chiral domain in the phase diagram shrinks 
as K is increased. In Fig. [H the first order phase tran- 
sition line is shown for three values of the Debye length 
for Nr — 100 and Nq = 3. It can be seen that even for 
= L the value of Uc barely changes as compared to its 
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K = value of 0.4 L. Our results suggest that ~ Oc 
sets the threshold where screening can wipe out the chiral 
ordering altogether. 

We have neglected the effect of counterions in the 
above study. This can be justified by noting that the 
polyelectrolytes are very short rods and unless we choose 
extremely high values for q, it is unlikely that they con- 
dense on the rods. This means that they will move rel- 
atively freely in the vicinity of the linear polyelectrolyte 
array, and merely participate in the Debye screening of 
the interaction. This effect can be taken into account by 
considering a local effective value for the Debye length 
that incorporates the density of the counterions as well 
as that of salt. Moreover, the presence of the counteri- 
ons in an effective "cell" around the array provides the 
necessary neutralization of the system. In our treatment, 
however, we did not need to care about the neutraliza- 
tion criterion because the centers of the charged rods 
are fixed in space and in effect it is only multipoles (not 
monopoles) that are interacting. 

Finally, we note that the results obtained here are 
related to the studies of the equilibrium configuration 
of ions in a confined plasma, where even similar chiral 
phases are observed as the confinement potential becomes 
anisotropic [l7j . 
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APPENDIX A: END-GRAFTED CHARGED 
RODS 

In this Appendix, we consider the case where the 
charged rods are actually grafted at their ends instead 
of the mid-points. This will make a model for a poly- 
electrolyte comb, i.e. a backbone with a linear array of 
charged side chains, which has in fact been synthesized 
and studied 

We have performed simulations on such a system and 
have found the same general behavior as the previous 
case. The corresponding phase diagram for the case of 
Nq = 7 is shown in Fig. [51 where the transition points 
are found to be Oc = 1.40 L, and oi = 0.3 L. Note that 
the region (b) in Fig. [5] corresponds to rods that are 
anti-parallel to their neighbors, which is a direct con- 



sequence of the non-vanishing dipole moment of rods 
that are grafted at one end as opposed to those that 
are grafted from the mid-point. Consequently, region (b) 
in this case corresponds to dipoles in a staggered field. 
Similarly, in the chiral phase the neighboring rods are 
out of phase by n (instead of tt/2) plus a residual twist 
angle. We also note that the transitions are occurring 



Disordered 



2.1 
1.3 



, Partially ordered 0=— 
pi} 2 



1A=0 (b) 




FIG. 9: Phase diagram of the polyelectrolyte comb model 
in the a ~ T plane. Dashed line: The crossover transition 
line from a disordered regime to partially ordered one [region 
(a)]. Dash-dotted line: The crossover transition line from 
a partially ordered regime to the staggered ground-state of 
the low-density regime [region (b)]. Solid line: First order 
transition line from a partially ordered regime to the chiral 
ordered phase of the intermediate-density case [region (c)]. 
Hashed region: The domain in which a multitude of transi- 
tions occur to phases with various kinds of orderings. Region 
(d): The "sea urchin" phase. Schematic configuration of rods 
corresponding to the different parts of the phase diagram is 
indicated in the lower panel. The values of the various tran- 
sition points indicated in the phase diagram are ac = 1.40 L 
and ai — 0.3 L. 



at larger values of the reduced temperature and lattice 
spacing, because of the stronger dipolar interactions. 
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